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THE  LIU,  KATSAROS,  AND  BUSINGER  (1979) 

BULK  ATMOSPHERIC  FLUX  COMPUTATIONAL  ITERATION 
PROGRAM  IN  FORTRAN  AND  BASIC 


ABSTRACT 

The  computer  program  described  by  Liu,  Katsaros,  and  Businger  (1979)  for  calculating  the  bulk- 
derived  atmospheric  fluxes,  stability,  and  roughness  is  presented  in  both  FORTRAN  and  BASIC  versions. 

The  results  of  12  test  calculations  are  presented  to  ensure  that  the  program  has  been  properly  encoded. 

INTRODUCTION 

Liu,  Katsaros,  and  Businger  (1979)  present  a  model  developed  for  the  marine  atmospheric  surface 
layer  which  takes  into  consideration  the  interfacial  sublayers  on  both  sides  of  the  air-sea  interface  where 
molecular  constraints  on  transports  are  important.  Flux-profile  relations  which  are  based  on  the  postu¬ 
lation  of  intermittent  renewal  of  the  surface  fluid  are  matched  to  the  logarithmic  profiles  and  compared 
with  both  field  and  laboratory  measurements.  These  relations  enable  numerical  determination  of  air-sea 
exchanges  of  momentum,  heat,  and  water  vapor  (or  bulk  transfer  coefficients)  employing  the  bulk 
parameters  of  mean  wind  speed,  temperature,  and  humidity  at  a  certain  height  in  the  atmospheric  sur¬ 
face  layer  and  the  water  temperature. 

With  increasing  wind  speed,  the  flow  goes  from  smooth  to  rough  and  the  bulk  transfer  coefficient 
for  momentum  also  increases.  The  increase  in  roughness  is  associated  with  increasing  wave  height, 
which  in  the  present  model  results  in  sheltering  at  the  wave  troughs.  Due  to  the  decrease  in  turbulent 
transport,  the  transfer  coefficients  of  heat  and  water  vapor  decrease  slightly  with  wind  speed  after  the 
wind  speed  exceeds  a  certain  value.  The  bulk  transfer  coefficients  are  also  found  to  decrease  with 
increasing  stability.  By  including  the  effects  of  stability  and  interface  conditions  in  bulk  parameteriza¬ 
tion,  the  model  provides  a  way  to  account  for  physical  conditions  which  are  known  to  affect  air-sea 
exchanges. 

A  computer  program  utilizing  computational  iteration  is  required  in  order  to  solve  three  simul¬ 
taneous  equations  of  three  unknowns,  which  are  formed  from  five  equations  containing  twelve  depen¬ 
dent  variables.  Although  the  original  program  was  written  in  FORTRAN,  the  recent  proliferation  of 
compact  desk-top  computers  usable  in  the  field  has  warranted  its  translation  into  BASIC. 

The  program  requires  seven  input  parameters: 

1.  Wind  speed  («)  in  m/s; 

2.  Air  temperature  ( 7)  in  °C; 

3.  Specific  humidity  (?)  in  kg/kg; 

4.  Water  temperature  at  the  surface  ( 7$)  in  °C; 

5.  Altitude  of  the  wind-speed  measurement  (zH)  in  m; 

6.  Altitude  of  the  air-temperature  measurement  ( z ,)  in  m;  and 

7.  Altitude  of  the  humidity  measurement  (z,)  in  m. 
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If  the  humidity  is  to  be  inputted  in  terms  of  the  wet  bulb  temperature,  dew  point  temperature,  or  rela¬ 
tive  humidity,  consult  the  Appendix  at  the  end  of  this  report. 

The  results  are  outputted  in  the  form  of  eight  parameters: 

1.  The  friction  velocity  (u.)  in  m/s; 

2.  The  scaling  potential  temperature  (0  ,)  in  °C; 

3.  The  scaling  specific  humidity  (?,)  in  kg/kg; 

4.  The  roughness  length  (z0)  in  m; 

5.  The  Monin-Obukhov  stability  (z/L); 

6.  The  wind-speed-roughness  Reynold’s  number  (Jtr); 

7.  The  temperature-roughness  Reynold’s  number  ( R ,);  and 

8.  The  humidity-roughness  Reynold’s  number  (Rq). 

From  the  outputted  information  it  is  possible  to  compute  the  momentum  flux  (A/)  in  N/m2,  the 
sensible  heat  flux  ( Hs )  in  W/m2,  and  the  humidity  flux  (£)  in  kg/  (s  m2)  by 

A/  —  —  pu\f 

HS-~P  Cp9t  u„ 

and 

E  -  -p  q,u¥ 

where  the  density  of  moist  air  (p)  in  kg/m3  and  the  specific  heat  of  moist  air  at  constant  pressure  (Cp) 
in  J/(kg/K)  can  be  determined  by 

(3.4838  x  1Q-3)F 
P  "  Tv  +  273.16 

and 

C,-  1.00411  +  0.9(g)]  x  103, 

where  P  is  the  barometeric  pressure  in  pascals  (if  not  recorded,  assumed  to  be  1.01325  x  10s  Pascals) 
and  Tv  is  the  virtual  potential  temperature  in  °C  such  that 

Tv  -  {(T  +  273.16)  x  [1  +  (0.608?)]}  -  273.16. 

By  convention,  a  positive  sign  is  used  to  indicate  an  upward  flux  and  a  negative  sign  a  downward  flux. 
Because  momentum  flux  is  almost  always  downward  (except  in  a  few  cases  when  the  water  velocity 
exceeds  the  wind),  the  negative  momentum  flux  is  traditionally  defined  as  the  stress  (r), 

r  =  -  M. 

Additionally,  it  is  frequently  convenient  from  a  thermodynamic  perspective  to  represent  the  humidity 
flux  in  terms  of  the  latent  heat  flux  ( HL )  in  W/m2  such  that 

HL  =  E  Ly 

where  L,  is  the  latent  heat  of  vaporization  in  J/kg, 

Ly  -  4.1868(597.31  -  0.56525  T)  x  103. 
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PROGRAM  LKB79F 

COMMON /PIN/U,T,Q,TS,ZU,ZT,ZQ, ID 
COMMON /POUT/USR  ,TSR , QSR , ZO ,ZL ,RR , RT , RQ 
C 

C  INPUT  DATA  FROM  TTY 
C 

1000  CONTINUE 
USR=0.0 
TSR=0 . 0 
QSR=0 . 0 
Z0=0.0 
ZL=0 . 0 
RR=0.0 
RT=0.0 
RQ=0.0 
CD=0.0 
WRITE  (5,1) 

1  FORMAT ( '  INPUT  VALUES  ARE: ' ) 

WRITE  (5,2) 

2  FORMATC  U=  ',  $) 

READ(5,900)  U 
WRITE(5,3) 

3  FORMATC  T=  ',  $) 

READ(5,900)T 

WRITE(5,4) 

4  FORMATC  Q=  ',$) 

READ(5,900)Q 

WRITE(5,5) 

5  FORMAT ( '  TS=  ’ ,  $ ) 

READ(5,900)TS 

WRITE(5,6) 

6  FORMATC  2U=  '  ,  $) 

READ(5,900)ZU 

WRITE(5,7) 

7  FORMAT ( '  ZT=  ' ,  $  ) 

READ(5,900)ZT 

WRITE(5,8) 

8  FORMATC  ZQ=  ',$) 

READ(5,900)ZQ 
CALL  ASL(IER) 

WRITE(5,20)USR,TSR,QSR,Z0,ZL,RR,RT,RQ,IER 
20  FORMATC  OUTPUT  VALUES  ARE:',/,'  USR=  '  .G13.6,/, 

A  '  TSR=',G13.6,/, 

B  '  QSR=',G13.6,/, 

C  '  ZO  =' ,G13.6,/, 

D  '  ZL  =' ,G13.6,/, 

E  '  RR  =',G13.6,/, 

F  '  RT  =*,G13.6,/, 

G  ’  RQ  =' ,G13.6,/,'  IER=',G13.6) 

900  FORMAT(G13.6) 

WRITE( 5,910) 

910  FORMATC  LAST  CASE?  0=YES,  l=NO  ',$) 

READ(5 ,920) IEND 
920  FORMAT(ll) 

IF  (IEND.EQ.I)GO  TO  1000 
END 

SUBROUTINE  ASL(IER) 

C 

C  TO  EVALUATE  SURFACE  FLUXES,  SURFACE  ROUGHNESS,  AND  STABILITY  OF 
C  THE  ATMOSPHERIC  SURFACE  LAYER  FROM  BULK  PARAMETERS  ACCORDING  TO 
C  LIU  EL  AL.  (T9)  JAS  36  1722-1735 
C  WRITTEN  BY  TIM  LIU  ON  5/8/79 
C 

C  INPUT: 

C  U  WIND  SPEED  IN  M/S 
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00660 
00670 
00680 
00690 
00700 
00710 
00720 
00730 
0071*0 
00750 
00760 
00770 
00780 
00790 
00800 
00810 
00820 
00630 
0081*0 
00850 
00860 
00870 
00880 
00890 
00900 
00910 
00 920 
00930 
009^0 
00950 
00960 
00970 
00980 
00990 
01000 
01010 
01020 
01030 
010U0 
01050 
01060 
01070 
01080 
01090 
01100 
OHIO 
01120 
01130 
011l*0 
01150 
01160 
01170 
01180 
01190 
01200 
01210 
01220 
01230 
012U0 
01250 
01260 
01270 
01280 
01290 
01300 
01310 


C  T  TEMPERATURE  IN  DEG  C 
C  Q  SPECIFIC  HUMIDITY  IN  KG/KG 
C  TS  SURFACE  TEMPERATURE  IN  DEG  C 
C  ZU  HEIGHT  OF  WIND  FENSOR  IN  METERS 
C  ZT  HEIGHT  OF  TEMPERATURE  SENSOR 
C  ZQ  HEIGHT  OF  HUMIDITY  SENSOR 

C  ID  SEE  SUBROUTINE  DRAG  FOR  DETAIL  DEFINITION ,ID=1  (KONDO), 
C  ID+2  ( SMITH), ID=3  (LARGE) 

C  OUTPUT: 

C  USR,TSR,QSR  SCALING  QUANTITIES  FOR  U,T,Q 
C  Z0,ZL  ROUGHNESS  AND  STABILITY  PARAMETERS 
C  RR,RT,RQ  ROUGNESS  REYNOLD  -NUMBERS  FOR  U,T,Q 
C 

C  IER=1  FAIL  TO  CONVERGE 
C  IER=2  LBK  ERROR 

C 

COMMON/PIH/U,T,Q,TS,ZU,ZT,ZQ,ID 
COMMON/ POUT/ USR , TSR ,QSR , ZO , ZL , RR , RT , RQ 
IER=0 

RI=9.8l»ZU*(T-TS ) / ( (273. 15+T ) *U**2 ) 

IF(RI.GT.0.25)IER=-1 

VISA=.15E-U 

ZL=0. 

Z0=.0005 

US=0. 

CALL  HUMLOW(TS,TS,QS) 

DU=U-US 

DT=T-TS 

DQ=Q-QS 

USR=.01*«DU 

N3=0 

30  CONTINUE 

N1=0 

10  CONTINUE 

U10=USR*AL0G ( 10 . /ZO )/.U 
TYPE  81*00, ID, U10, CD 
81*00  FORMAT  ('  DRAG  ’  .3G13.6) 

CALL  DRAG(ID,U10,CD) 

TYPE  81*00,  ID, U10, CD 
C=1./SQRT(CD) 

Z0N=10./EXP(  .1**C) 

TEST1=ABS ( ( ZON-ZO ) / ( ZO+1 . E-8 ) ) 

TYPE  88lO,TESTl ,ZON,ZO,C ,CD,N1 
8810  FORMAT  ('  TEST1 ,ZON,ZO,C ,CD,N1 ' ,6G13.6) 
IF(TESTl.LT.O.Ol)GO  TO  19 
N1=N1+1 

IF(N1.GT.50)G0  TO  95 
ZO=ZON 
GO  TO  10 
19  CONTINUE 

PUZ=PSI ( 1 ,ZL) 

ZTL=ZL*ZT/ZU 
ZQL=ZL*ZQ/ZU 
PTZ=PSI ( 2 , ZTL ) 

PQZ=PSI(2,ZQL) 

USR=DU*0.1*/(  ALOG  ( ZU/  ZO )  -PUZ ) 

RR=ZO*USR/VISA 
ZTSR=ZT*USfl/VISA 
ZQSR=ZQ*USR/VISA 
CALL  LKB(RR,RT,1) 

IF(RT  .NE.  -999* )G0  TO  21 
IER=2 

WRITE(5,2)RR 

RETURN 

21  CALL  LKB(RR,RQ,2) 

IF( RT  . NE.  -999. )G0  TO  22 
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013  00 

IER-2 

01330 

WRITK(5,2)RR 

013*40 

RETURN 

01350 

22 

S=2 . 2  * ( ALOG ( ZTSR / RT ) -PTZ ) 

01360 

P=2.2*(AL0G(ZQSR/RQ)-PQZ) 

01370 

TSR=DT/S 

01380 

QSR=DQ/D 

0.1390 

CALL  Z  ETA ( T , Q , US  R ,TSR ,QSR , ZU , ZLN ) 

01*100 

TEST3=ABS ( ( ZL-ZLN ) / ( ZL+1 . E-8 ) ) 

01*410 

IF (TEST3.LT. 0.01)00  TO  39 

01*120 

N3=N3+1 

01*430 

IF(N3.GT. 50 )G0  TO  95 

01*4*40 

ZL=ZLN 

01*4  50 

GO  TO  30 

Oi*i6o 

39 

CONTINUE 

01*470 

GO  TO  90 

01*430 

95 

IER=1 

01  *4  90 

WRITE(5,1)N1,N3 

01500 

1 

F0RMAT(1X,21HASL  FAILS  TO  CONVERGE, 315) 

01510 

2 

FORMAT( IX ,21HLKB  FAILS  BECAUSE  RR=,E12.U) 

01520 

99 

RETURN 

01530 

END 

015*40 

SUBROUTINE  DRAG( ID,U,CD) 

01550 

C 

01560 

C  TO  DETERMINE  NEUTRAL  DRAG  COEFFICIENT  CD  FROM  WIND  SPEED 

01570 

C  AT  1C 

1  M  U  IN  M/S 

01580 

C  ID=1 

K0ND0(1975)  BLM  9  91-112 

01590 

C  ID=2 

SMITH ( 1980)  JPO  10  709-726 

01600 

C  ID=3 

LARGE  &  POND  (1981)  JPO  11  324-336 

01610 

C  RANGE 

:  OF  U  SPECIFIED  ARE:  KONDOt . 3 ,50 ) ,SMITH(6 ,22 ) ,L&P( 4 ,25 

01620 

C  WRITTEN  BY  TIM  LIU  FOR  VAX  ON  2/10/82 

01630 

C 

016*40 

DIMENSION  RAN(5) ,A(5) ,B(5) ,P(5) 

01650 

DATA  RAN/2.2,5.,8.,25.,50./ 

01660 

DATA  A/0,,0.771,0.867,1.2,0./ 

01670 

DATA  B/ 1.08, 0.0858,. 0667, 0.025, 0.073/ 

01680 

DATA  P/-0.15.1. ,1. ,1. ,1./ 

01690 

K=ID-2 

01700 

IF('K)100,200,300 

01710 

100 

IF(U .GT . 50 . )G0  TO  131 

01720 

IFtU.LT. -3)G0  TO  130 

01730 

1=1 

017*40 

110 

IF(U. LE.  RAN ( I ) )G0  TO  120 

01750 

1  =  1  +  1 

01760 

GO  TO  110 

01770 

120 

CD=(A(I)+B(I)*U»*P(I) )/1000. 

01780 

GO  TO  99 

01790 

130 

CD=1 . 5E-03 

01800 

GO  TO  99 

01810 

131 

CD=3 .7E-03 

01820 

GO  TO  99 

01830 

200 

CD=(0.6l+0.063*U)/1000. 

olOlo 

GO  TO  99 

01850 

300 

IF(U.LT.ll. )G0  TO  301 

01860 

CD=(0 . 49+0,o65*U )/1000. 

01870 

GO  TO  99 

01880 

301 

CD=1 . 2E-3 

01890 

99 

RETURN 

01900 

END 

01910 

SUBROUTINE  HUMLOW(T,TW,Q) 

01920 

c 

01930 

C  TO  EVALUATE  SPECIFIC  HUMIDITY  Q  FROM  DRY  AND  WET  BULB  TEMP 

019*40 

C  T  AND  TO  ACCORDING  TO  LOWE(77)  JAM  16  100-103 

01950 

C  WRITTEN  BY  TIM  LIU  ON  5/3/79,  REVISED  FOR  VAX  ON  2/10/82 

C1960 

C 

01970 

DIMENSION  A( 6) 

5 


01960 

DATA  A/4.4 365I9E-1 , 1 . 428946E-2 ,2 . 650649E-4 , 3. 031240E-6 , 

0 1990 

$  2.03408lE-8,6. 136821E-11/ 

02000 

P=1013.25 

02010 

X=0. 

02020 

DO  100  1=1,6 

02030 

J=7-I 

02040 

X=(X+A( J ) )*TW 

02050 

100 

CONTINUE 

02060 

ES=6.107800+X 

02070 

W=0 . 622»ES/ ( P-ES ) -4 , 045E-04*(T-TW ) 

070(30 

RETURN 

07090 

END 

02100 

SUBROUTINE  LKB( RR,RT, IFLAG ) 

U,’  1  I  'J 
0J1J0 

c 

0  TO 

DETERMINE  THE  LOWER  BOUNDRY  VALUE  RT  OF  THE  LOGARITHMIC 

OJl  iU 

C  PROFILES  OF  TEMPERATURE  ( IFLAG=1 )  OR  HUMIDITY  (IFLAG=2) 

Oi'lU  0 

C  IN 

THE  ATMOSPHERE  FROM  ROUGHNESS  REYNOLD  NUMBER  RR  BETVEEN 

OJ  l‘»0 

C  0  AND  1000.  OUT  OF  RANGE  RR  INDICATED  BY  RT=-999. 

0  J.  iou 

0  BASED  ON  LIU  ET  AL. .  (1979)  JAS  36  1722-1723 

OJl'fU 

C  WRITTEN  BY  TIM  LIU  ON  1/22/78,  REVISED  FOR  VAX  ON  2/10/62 

OJltfG 

C 

uji;o 

DIMENSION  A(8,2),B(8,2) -RAN (8 ) 

02. '00 

DATA  A/0. 177, 1.376, 1.026, 1.625, 4. 66l, 34, 904, 1667. 19, 5, 88E5, 

02710 

$  O.292 , 1.808,1. 393, 1-956, 4. 991*  ,30.709,  l^S.  68,2. 98E5/ 

02220 

DATA  B/0. ,0.929,-0.599,-1.018,-1.475,-2.067.-2.907,-3.935, 

0  2230 

$  0. ,0.826,-0.528,-0.870,-1.297,-1.845,-2.682,-3.616/ 

02240 

DATA  RAN/O. 11, 0.825, 3. 0,10. 0,30. 0,100. ,300. ,1000./ 

02750 

1  =  1 

02260 

IF  ( RR. LE.O. .OR.RR.GE. 1000. )  GO  TO  90 

02270 

10 

CONTINUE 

07260 

IF  (RR.LE.RAN(I) )  GO  TO  20 

02290 

1=1+1 

02  300 

GO  TO  10 

02310 

20 

RT=A( I , IFLAG ) »RR»»B( I , IFLAG ) 

07320 

GO  TO  99 

02330 

90 

HT=-999. 

02340 

99 

RETURN 

07330 

END 

02  360 

FUNCTION  PSI(IU,ZL) 

02370 

C  TO 

EVALUATE  THE  STABILITY  FUNCTION  PSI  FOR  WIND  SPEED  (IFLAG=1) 

02  360 

C  OR 

FOR  TEMPERATURE  AND  HUMIDITY  PROFILES  FROM  STABILITY  PARAMETER 

02390 

0  SEE  LIU  ET  AL.  (1979)  JAS  36  1722-1723  FOR  DETAILS 

02400 

C  WRITTEN  BY  TIM  LIU  ON  9/12/71,  REVISED  FOR  VAX  ON  2/10/82 

02410 

C 

02420 

IF( ZL) 10 ,20 , 30 

02430 

10 

CHI=(l.-l6.*ZL)**0,25 

02440 

I F ( ID.EQ. 1 )G0  TO  11 

02450 

PSI=2.*AL0G( (l.+CHI*CHl)/2.) 

02460 

GO  TO  99 

02470 

11 

PSI=2.*AL0G( (l.+CHl)/2. )+ALOG( ( 1 , +CHI*CHI ) /2. )-2.*ATAN(CHI ) 

02430 

&  +2 . *ATAN ( 1 . ) 

02490 

GO  TO  99 

02500 

20 

PS  1=0. 

02510 

GO  TO  99 

02570 

30 

PSI=-6.*AL0G(1.+ZL) 

073  30 

99 

RETURN 

02340 

END 

07350 

SUBROUTINE  ZETA(T,Q,USR,TSR,«SR,Z,ZL) 

023o0 

C 

07370 

C  TO 

EVALUATE  OBUKHOVS  STABILITY  PARAMETER  Z/L  FROM  AVERAGE 

02530  C  TEMP  T  IN  DEG  C,  AVERAGE  HUMIDITY  Q  IN  GM/GM,  HEIGHT  Z  IN  M, 

02390  C  AND  FRICTIONAL  VEL.TEMP.  ,  HUM.  INMKS  UNITS 

u:’6uO  C  GEE  LIU  ET  AL.  (1979)  JAS  36  1722-1723  FOR  DETAILS 

021.10  C  WRITTEN  BY  TIM  LIU  ON  10/1/77,  REVISED  FOR  VAX  ON  2/10/62 

02020  C 

020.50  VONO.L 
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0261*0  G=9.8l 

020%  TA=2','0.  lo+T 

CVooC  TV=TA»(l.+0.6l*Q) 

02670  TVSR=TSR*  tl . +0 . 6l*Q ) +0 . 6l«TA»QSR 

Oot.rtO  IF(TVSR.EQ.0.)G0  TO  10 

0260 J  0B=TV*USR*U3R/ (G»V0N*TVSR ) 

D/XO  2L=Z/0B 

02  '  *  0  GO  TO  99 

0  j"7.'0  1C  ZL=0. 

jj'Mj  ,)  RETURN 

Oc'VU'j  END 


THE  BASIC  PROGRAM 


10 
20 
30 
40 
50 
60 
70 
80 
90 
100 
1  10 
1  20 
130 
140 
150 
160 

1  70 
180 
190 
200 
210 
220 
230 
240 
250 
260 

2  70 
280 


nen 
REM 
REM 
Ren 
Ren 
Ren 
Ren 
Ren 
Ren 
Ren 
REM 
REM 
REM 
REM 
REM 
REM 
REM 
REM 
REM 
REM 
DIM  R 
DATA 
DATA 
DATA 
DATA 
DIM  A 
DATA 
DATA 


"LKB79B"  PROGRAM 


12  OCT  83 


T  V  BLANC 


NRL  4110 


DISK  #6 


THE  BULK  FLUX  COMPUTATIONAL  ITERATION  PROGRAM  OF  LIU,  KATSAROS,  t- 
BUS I NCER  (1979)  TRANSLATED  FROM  FORTRAN  INTO  ELEMENTRY  BASIC  FOR  USE 
ON  A  HEWLETT-PACKARD  MODEL  9845  THIS  PROGRAM  WAS  ADAPTED  FROM  AN 
0RG1NAL  TRANSLATION  OF  A  T  WILSON  NRL  2820  FOR  USE  ON  A  DIGITAL 
EQUIPMENT  CORPORATION  MODEL  DEC-10  COMPUTER 


INPUTS 


NOTES 


U=-  Wind  Speed  (m/s),  T-  Air  Temp  <C>,  Q=  Spec  Humid  (Kg/Kg). 
Tl=Tui=  Water  Temp  <C),  21=2u=  U  Altitude(m),  Z2cZt=  T  Alt  (m) 

Z3=Zq=  Q  Alt  <m> 


U  SHOULD  NOT 


0 


OUTPUTS  U 1 -Us  tr  =  U*  (m/s),  T2=Tstr=  T*  (C>.  Ql=Qstr=  Q*  (Kg/Kg), 

Z4=7o-  Rough  Length  (m).  Z6=Z/L=  Mon l n-Ob u k h o v  Stability, 

Rl=Ru-  U  Rough  Reynolds  No  ,  R2=R  t  =  T  Rough  Reynolds  No 
R3=Rq=  Q  Rough  Reynolds  No 

5(5),  A2 ( 5 ) ,  DO ( 5 ) ,  P ( 5  ) 

2  2.  5  ,  8  ,  25  ,  50 

0  ,0  7/1,0  86  7,  1  2,  O 
1  08.0  0858,0  0667,0  025,0  073 
-0  15,  1  ,1  ,1  ,1 

1(6) 

4  4365 1 9E -1,1  428946E-2.2  650649E-4 

3  031240E-6,  2  034081E-B. 6  136821E-11 


290  DIM  A<8,  2),  B(B.  2),  R6 ( 8 ) 

300  DATA  0  177,1  376,1  026.1  625.4  661,34  904,1667  19,5  88E5 

310  DATA  0  292,  1  BOB.  1  393.  1  956,4  994,30  709,  1448  68,2  98E5 

320  DATA  0  .O  929,-0  599,-1  01B.-1  475,-2  067,-2  907,-3  935 

330  DATA  O  .0  826,-0  528,-0  870,-1  297,1  845,-2  682,-3  616 

340  DATA  0  11,0  825,3  0,10  0,30  0,100  0,300  ,  1 OOO  O 
350  RESTORE  220 
360  REM 

370  FOR  1=1  TO  5 
380  READ  R5(I) 

390  NEXT  I 

400  FOR  1=1  TO  5 

410  READ  A2( I ) 

420  NEXT  I 

430  FOR  1=1  TO  5 

440  READ  B2 ( I ) 

450  NEXT  I 

460  FOR  1=1  TO  5 

470  READ  P ( I ) 

480  NEXT  I 

490  FOR  1=1  TO  6 

500  READ  Al(I) 

510  NEXT  I 
520  FOR  J=1  TO  2 
530  FOR  1=1  TO  8 
540  READ  A( I,  U) 

550  NEXT  I 
560  NEXT  J 
570  FOR  J=1  TO  2 
580  FOR  1=1  TO  8 
590  READ  B< I,  J) 

600  NEXT  I 
610  NEXT  J 
620  FOR  1=1  TO  8 
630  READ  R6( 1 ) 

640  NEXT  I 
650  REM 
660  PRINT  " 
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670  PRINT  " 

680  PRINT  "INPUT  DATA:  U ( m/ s  > ,  T  <  C  ) ,  0 <  k  g  / k  g  ) ,  Tin <  C  > ,  Z u  <  m  > ,  Z  t  <  m ) ,  Zq  <  m ) 

690  LET  U1=0 
700  LET  T2=0 
710  LET  01*0 
720  LET  Z4=0 
730  LET  Z6»0 
740  LET  R 1 B0 
750  LET  R2-0 
760  LET  R3<=0 

770  INPUT  U, T, Q, Tl, Zl, Z2, Z3 
780  PRINT  U, T, 0, Tl, Zl, Z2, Z3 

790  REM  .  ....  CALL  SUROUTINE  ASL(IER) 

800  LET  A9=I9 
810  COSUB  990 
820  LET  I9=A9 
830  REM 


840 

PRINT 

"OUTPUT  VALUES  ARE 

850 

PRINT 

"U*  <m/*> 

=  » 

,  U1 

860 

PRINT 

"T*  (C) 

ac  " 

,  T2 

870 

PRINT 

"0*< kg/kg ) 

as  " 

.  01 

880 

PRINT 

"Zo  (m) 

,  Z4 

890 

PRINT 

"Z/L 

,  Z6 

900 

PRINT 

"Ru 

*  " 

,  R 1 

910 

PRINT 

"Rt 

,  R2 

920 

PRINT 

"Rq 

sc  " 

,  R3 

930  PRINT  " 

940  REM 

950  PRINT  "LAST  CASE^  0=YES,  1=N0  "; 

960  INPUT  18 

970  IF  18=1  THEN  GOTO  660 
980  GOTO  2670 

990  REM  . SUBROUTINE  ASL(IER) 

1000  LET  A9=0 

1010  LET  R4=9  61#Z1*<T-T1 >/< (273  1 5+T ) *U'2 ) 

1020  IF  R4>.  25  THEN  GOTO  1040 

1030  GOTO  1050 

1040  LET  A9=- 1 

1050  LET  Vl=l  5E-5 

1060  LET  Z5=0 

1070  LET  Z4=  0005 

1080  LET  U2=0 

1090  LET  H1=T1 

1100  LET  H2=T 1 

1 1 10  LET  H3=01 

1120  REM  CALL  HUMLOW <  Tin,  Tin,  OS  > 

1130  GOSUB  2290 
1140  LET  T 1 =H1 
1150  LET  T 1 =H2 
1 160  LET  Q1=H3 
1170  LET  D0=U-U2 
1180  LET  WP=T-T1 
1190  LET  W3=Q-01 
1200  LET  U1  *=  04*D0 
1210  LET  N3=0 

1220  REM  CONTINUE 

1230  LET  N1=0 

1240  LET  U0=U1*L0G( 10/Z4)/  4 
1250  LET  D1=D1 
1260  LEI  D2=U0 
1270  LET  D3=C 1 
1280  GOSUB  2050 
1290  REM 
1300  LET  D1=D1 
1310  LET  UO=D2 
1320  LET  C1=D3 


CALL  DR AG ( ID. U10, CD  > 


1330 

1340 

1350 

1360 

1370 

1380 

1390 

1400 

1410 

1420 

1430 

1440 

1450 

1460 

1470 

1480 

1490 

1500 

1510 

1520 

1530 

1540 

1550 

1560 

1570 

1580 

1590 

1600 

1610 

1620 

1630 

1640 

1650 

1660 

1670 

1680 

1690 

1700 

1710 

1720 

1730 

1740 

1750 

1760 

1770 

1780 

1790 

1800 

1810 

1820 

1830 

1840 

1850 

I860 

1870 

1880 

1890 

1900 

1910 

1920 

1930 

1940 

1950 

1960 

1970 

1980 


LET  C=1/SQR(C1) 

LET  Z5=10/EXP<  4*0 

LET  T3=ABS< < Z5-Z4 ) / < Z4+1E-B) ) 

IF  T3C  01  THEN  GOTO  1410 
LET  N1=N1+1 

IF  Nl>50  THEN  GOTO  2020 
LET  Z4=Z5 
GOTO  1240 

REM  .  CONTINUE 

LET  P1=FNP< 1. Z6> 

LET  Z7=Z6*Z2/Z1 
LET  Z8=Z6*Z3/Z 1 
LET  P2=FNP (2.  Z7) 

LET  P3=FNP ( 2.  Z8> 

U1=D0*  4/ <  LOG  ( Z 1 /Z4 )-P 1 ) 

LET  R1=Z4*U1/V1 
LET  Z9=Z2*U1/V1 
LET  Z0=Z3*U1/V1 
LET  L1=R1 
LET  L2=R2 
LET  L3=*l 
GOSUB  2410 

REM  . .  CALL  LKB  (Ru.  fit.  1  ) 

LET  R 1 =L 1 
LET  R2=L2 
LET  L3=L3 

IF  R20-999  THEN  GOTO  1630 
LET  A9=2 

PRINT  "LBK  FAILS  BECAUSE  Ru=",Rl 
RETURN 

REM  SET  VALUE  OF  2  TO  A  VARIABLE 

LET  L1=R1 

LET  L2=R3 

LET  L3=2 

GOSUB  2410 

LET  R1=L1 

LET  R3=L2 

LET  L3=L3 

IF  R20-.  999  THEN  GOTO  1740 
LET  A9=2 
GOTO  1610 

LET  S=2  2*  <  LOG  <  Z9/R2 ) -P2 ) 

LET  D=2.  2*<L0G< Z0/R3J-P3) 

LET  T2=W2/S 
LET  Q1=W3/D 

REM . CALL  ZETAtT, Q. Ustr,  Tstr, Qstr,  ZU,  ZLN) 

LET  F1=T 
LET  F2=Q 
LET  F3=U1 
LET  F4=T2 
LEI  F5-Q1 
LET  F6=Z1 
LET  F7= Y 1 
GOSUB  2550 
LET  T=F1 
LET  Q=F2 
LET  U1=F3 
LET  T2=F4 
LET  G1=F5 
LET  Z 1=F6 
LET  Y 1 =F7 

LET  T4=*ABS(  <Z6-Y1  )/<Z6+lE-8)  ) 

IF  T4<.  01  THEN  GOTO  2000 
LET  N3=N3+1 

IF  N3>50  THEN  GOTO  2020 
LET  Z6=Y 1 
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1990  GOTO  1220 

2000  REM . CONTINUE 

2010  GOTO  2040 
2020  LET  A9= 1 

2030  PRINT  "ASL  FAILS  TO  CONVERGE" .  N1 ,  N3 
2040  RETURN 


2050  REM . SUBROUTINE  DRAG ( D1 .  D2.  D3) 

2060  REM  R5=RAN,  A2=A.  B2=B  IN  DRAG 


2070  LET  K=D 1 —2 
2080  IF  K— 0  THEN  GOTO  2220 
2090  IF  K>0  THEN  GOTO  2240 
2100  IF  D2>50  THEN  GOTO  2200 
2110  IF  D2<  3  THEN  GOTO  2180 
2120  LET  1=1 

2130  IF  D2<=R5< I )  THEN  GOTO  2160 
2140  LET  1=1+1 
2150  GOTO  2130 

2160  LET  D3=(A2(  I  >+B2<  I  )*D2''P<  I )  1/1000 

2170  GOTO  2280 

2180  LET  D3=l.  5E-3 

2190  GOTO  2280 

2200  LET  D3=3.  7E-3 

2210  GOTO  2280 

2220  LET  D3=C.  61+.  063*D2)/1000 
2230  GOTO  2280 

2240  IF  D2<1 1  THEN  GOTO  2270 
2250  LET  D3=(  .  49+. 065*D2)/1000 
2260  GOTO  2280 
2270  LET  D3=l.  2E-3 
2280  RETURN 


2290  REM . SUBROUTINE  HUMLOWIT.  TW,  Q) 

2300  REM  SUBROUTINE  HUMLOWtHl,  H2, H3) 

2310  REM  A1=A  IN  HUMLOW 


2320  LET  P=1013.  25 
2330  LET  X=0 
2340  FOR  1=1  TO  6 
2350  LET  J=7-I 
2360  LET  X=  <  X+Al <  J )  ) *H2 
2370  NEXT  I 

2380  LET  El=6  107800+X 

2390  LET  H3=  622*E1/(P-E1 )-4.  045E-4*<H1-H2) 

2400  RETURN 

2410  REM  . SUBROUTINE  LKB(Ru/Rt» IFLAG) 

2420  REM  SUBROUTINE  LKB < LI ,  L2.  L3 > 

2430  REM  R6=RAN  IN  LKB 

2440  LET  1=1 

2450  IF  L1<=0  THEN  GOTO  2530 
2460  IF  LI >=1000  THEN  GOTO  2530 

2470  REM . CONTINUE 

2480  IF  L1<=R6< I )  THEN  GOTO  2510 
2490  LET  1=1+1 
2500  GOTO  2470 

2510  LET  L2=A(I.L3)*R1'B(I,L3) 

2520  GOTO  2540 
2530  LET  L2=-999 
2540  RETURN 

2550  REM . . SUBROUTINE  ZETA  (  T.  Q.  Us tr ,  Tstr,  Ostr .  Z.  Z/L  > 

2560  LET  V=  4 
2570  LET  G=9  81 
2580  LET  T 5=273  16+F1 
2590  LET  T6=T5*(1+  61*F2> 

2600  LET  T7=F4*< 1+  61*F2)+.  61*T5*F5 
2610  IF  T7=0  THEN  GOTO  2650 
2620  LET  01=T6*F3*F3/(G*V*T7) 

2630  LET  F7=F6/01 
2640  GOTO  2660 
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2650  LET  F7=0 
2660  RETURN 
2670  END 

2680  REN  . 

2690  REN  .  FUNCTION  PSI  ( ID,  Z/L) 

2700  DEF  FNP( 17. X5> 

2710  REN . TEST  Z/L  FOR  0.  + 

2720  IF  X5=0  THEN  GOTO  2800 
2730  IF  X5>0  THEN  GOTO  2820 
2740  LET  C2=( 1-16*X5)~. 25 
2750  IF  17=1  THEN  GOTO  2780 
2760  RETURN  2*L0G ( ( 1 +C2*C2 ) /2 ) 

2770  GOTO  2830 

2780  RETURN  2*L0G<  < 1 +C2 ) /2  M-LOG ( < 1 +C2*C2 ) /2 > ~2*ATN< C2 > +2*ATN ( 1 > 

2790  GOTO  2830 

2800  RETURN  0 

2810  GOTO  2830 

2820  RETURN  -6*L0G ( 1 +X5 ) 

2830  FNEND 


TEST  CALCULATION  RESULTS 


To  ensure  that  the  programs  have  been  properly  copied  by  the  user,  the  results  of  12  test  calcula¬ 
tions  are  presented.  The  outputted  results  should  agree  with  those  indicated  below  to  within  ±2  of  the 
fourth  significant  figure. 


Test 

Input 

Output 

1 

u  =  9 

u„  -  0.3490 

-s 

1! 

K> 

0.  -  -6.922  x  10" 2 

q  =  0.007 

qm  -  -1.062  x  10~4 

Ts  =  14 

20  -  2.501  x  10“4 

■=N 

II 

o 

z/L  -  -9.937  x  10“2 

2,  —  10 

Rr  -  5.820 

o 

II 

R,  -  0.2705 
-  0.4226 

Test 

Input 

Output 

2 

u  =  9 

w,  -  0.3262 

T  =  14 

0,-  6.536  x  10~2 

q  =  0.007 

qt  -  -5.817  x  10~5 

Ts  -  12 

2Q  -  2.429  x  10~4 

o 

1 

z/L  -  7.087  x  10“2 

2,-  10 

Rr  -  5.282 

o 

II 

*4* 

R,  -  0.2986 

Rq  -  0.4598 

Test 

Input 

Output 

3 

</»  9 

u,  -  0.3514 

T  -  12 

0,- -6.981  x  10~2 

q  =  0.002 

qm  -  -2.879  x  10~4 

Ts  —  14 

20  -  2.501  x  10"4 

2U-10 

z/L  -  -0.1347 

2,-  10 

Rr  -  5.859 

z,-  10 

R,  -  0.2686 

R,  -  0.4201 

Test 

Input 

Output 

4 

u  -  9 

«,  -  0.3334 

7-  14 

0,  -  6.614  x  10"2 

q  -  0.002 

q ,  -  -2.297  x  10"4 

rs-  12 

20  -  2.466  x  10“4 

2„-  10 

z/L  -  3.175  x  10"2 

2,-  10 

Rr  -  5.481 

<4* 

1 

o 

R,  -  0.2875 

Rq  -  0.4452 

Test 

Input 

Output 

5 

u  -  2 

w,  -  7.230  x  10"2 

T-  12 

0,-  -9.516  x  10"2 

q  -  0.007 

q .  -  -1.486  x  10“4 

T,  -  14 

20  -  2.653  x  10"5 

2„-10 

z/L  -  -3.175 

2,  —  10 

Rr  -  0.1279 

2„-  10 

R,  -  0.2036 

Rq  -  0.3306 
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Test 

Input 

Output 

6 

u-2 

u,  -  3.731  x  10“2 

T-  14 

9,  -  4.336  x  10~2 

q  -  0.007 

q ,  -  -3.831  x  10~5 

7>  12 

z0  -  4.498  x  10-5 

zu  -  10 

zIL  -  3.578 

z,-  10 

R,  -  0.1119 

o 

1 

*,-0.1799 
*?  -  0.2962 

Test 

Input 

Output 

7 

u  -  18 

u.  -  0.7361 

T-  12 

fl,- -5.599  x  lO"2 

g  -  0.007 

q%  -  -8.643  x  10"5 

7;  -  14 

z0  -  5.289  x  10~4 

zu-  10 

z/L  -  -1.802  x  10"2 

z,-  10 

Rr  -  25.96 

z^-10 

R,  -  3.824  x  10"2 

R„  -  7.315  x  lO"2 

Test 

Input 

Output 

8 

u-  18 

u.  -  0.7259 

T-  14 

9  ,  -  5.542  x  10"2 

q  -  0.007 

-  -4.975  x  1(T5 

r,-  12 

z0  -  5.289  x  10~4 

zu-  10 

z/L  -  1.203  x  l(T2 

z,  -  10 

Rr  -  25.59 

O 

1 

*,  -  3.904  x  10“2 

Rq  -  7.449  x  10"2 

Test 

Input 

Output 

9 

«-  9 

u%  -  0.3250 

T-  12 

9,-  -7.050  x  lO'2 

q  -  0.007 

-  -1.081  x  lO'4 

rs  -  14 

z0  -  2.421  x  10"4 

O 

r^i 

1 

•4* 

z/L  -  -0.3487 

z,  -  10 

Rr  -  5.246 

1 

o 

*,  -  0.3007 

Rq  -  0.4625 

Test 

Input 

Output 

10 

w  -  9 

u,  -  0.3495 

r-  12 

0,  -  -6.613  x  10"2 

$  -  0.007 

q  -  -1.059  x  lO"4 

7;-  14 

z0  -  2.556  x  10"4 

zu  -  10 

z/L  -  -9.575  x  lO"2 

o 

CO 

1 

1  n" 

Rr  -  5.956 

o 

1 

►r 

R,  -  0.2642 
-  0.4142 
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Test 

Input 

Output 

11 

u  -  9 

u,  -  0.3489 

f-  12 

0,- -6.919  x  10“2 

q  -  0.007 

q ,  -  -1.016  x  10"4 

r,-  14 

20  -  2.501  x  10~4 

zu  —  10 

z/L  -  -9.777  x  10"2 

2,  -  10 

R,  -  5.818 

2,  -  30 

R,  -  0.2706 

R,  -  0.4227 

Test 

Input 

Output 

12 

u  —  6 

u,  -  0.2359 

T-l 

0.  -  -3.787  x  10"2 

q  -  0.004 

q,  -  -1.016  x  10“4 

T.  -  8 

20  —  1.536  x  10~4 

2U  -  5 

z/L  -  -7.009  x  10~2 

2,  —  15 

Rr  -  2.416 

2,-25 

R,  -  0.6049 

R,  -  0.8743 
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APPENDIX 


The  specific  humidity  ( q )  in  kg/ kg  is  calculated  by 

0.622e 

q  “  P  -  (0.378<?) 

where  P  is  the  barometric  pressure  and  e  is  the  water  vapor  pressure  in  pascals.  If  the  humidity  is 
measured  by  the  wet  and  dry  bulb  temperatures  (Twb,T)  in  °C  then 

e  -  e,  —  6.6  x  10"4[l  +  (1.15  x  10-3rwfc)]/>(7’ -  Twb) 

where  es  is  the  saturated  vapor  pressure  in  pascals  such  that 

es- Pat3  x  10-A  +  «4*5  +  -s*. 

where 

373.16 

”  r  +  273.16  ’ 
a2~  a i  -  1, 


fl4-  (10“j6>)-  1. 
a5  -  (10°3*2)  -  1, 

and  the  GofT-Gratch  humidity  formulation  constants  are 

bi -  3.49149,  -  7.90298, 

b2  -  11.344,  bs-  8.1328  x  10~3, 

*3  -  5.02808,  b6 - 1.3816  x  10"7. 

If  the  humidity  is  measured  by  the  dew  point  temperature  (7^)  in  °C  then 

e  -  Pat 3  x  10°2*4  +  °4*5  +  °s*6 

where 

373.16 
7^  +  273.16 

and  a2  through  a5  are  calculated  in  the  same  manner  as  for  the  saturated  vapor  pressure.  If  the  relative 
humidity  (RH)  in  %  is  measured  then 

e,RH 
6  ”  100 

where  the  saturated  vapor  pressure  ( es )  is  calculated  from  the  air  temperature  as  shown  above. 
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